######################################################################
# Time-stamp: <05/08/09 17:02:30 ostolop>
######################################################################
package EP::EPC::Core::Sequence_Motif;

use strict;

use EP::Config;
use EP::Common;
use EP::Common::General;
use EP::Common::Error qw ( :try );

use EP::EPC::Core;
use EP::EPC::Core::Component;
use EP::Core::Data::Parser::seq;

use EP::Core::Dataset::Data::seq;

use CGI qw/:standard/;
use File::Temp;
use File::Basename;

use Bio::DB::GFF;

use base qw ( EP::EPC::Core::Component );

sub process_component {
  my $self = shift;

  my $program = $self->{query}->param( "program" );
  my $result;

  if ( $program eq 'spexs' ) {
    $result = $self->_run_SPEXS ();
  } else {
    die "Do not know how to run this program: $program";
  }

  my $process_result = { spexs => $result,
			 src_dataset_id => $self->{src_dataset}->id,
			 dst_dataset_id => $self->{dst_dataset}->id
		       };	# HACK: produce HTML output and place it all directly in <html></html> as a result!

  $self->{CD}->{process_result} = $process_result;
  return $self->record_analysis_history ();
}

sub get_fg_dataset {
  my $self = shift;
  my @ids;

  my $src_dataset = $self->{src_dataset};
  my $selected_cluster = $self->{query}->param('current_cluster_selection');

  if ( $selected_cluster ) {
      # a cluster has been selected, get those ids
      @ids = $self->get_gene_ids_from_cluster( $selected_cluster );
      throw EP::Common::Error( -text => "No genes found for selected cluster" ) if scalar(@ids) == 0;
  } else {
      # get the genes from the currently selected dataset
      return $src_dataset if $src_dataset->type eq "seq";

      # Get the gene names ... for now just assume the NAMES entry in the row annotation is ok...
      @ids = map { $src_dataset->data->rows->{$_}{full_annot}->Values(0) }  sort { $a <=> $b } keys %{ $src_dataset->data->rows };      
      if ( scalar(@ids) == 0 || join("", @ids) =~ /^\s*$/ ) {
	  throw EP::Common::Error( -text => "No genes found for dataset nr $src_dataset" );
      }
  }

  # get the sequences for the cluster
  $self->{session}->status ( "Retrieving cluster sequences ..." );
  print STDERR "Going to retrieve ids: " . join( ", ", @ids ) . "\n";
  my $seq = EP::Core::Data::Parser::seq->_get_sequences ( $src_dataset->species, \@ids, $self->{query}->param('upstream_bp') );
  my $fg_dataset = $src_dataset->copy ( { name => undef } );
  $fg_dataset->type ( 'seq' );

  my $data = new EP::Core::Dataset::Data::seq ( sequences => $seq );

  $fg_dataset->data ( $data );
  $fg_dataset->data->start_bp( -1 * $self->{query}->param('upstream_bp') );
  $fg_dataset->data->end_bp( 0 );
  
  EP::Core::Data::Parser::seq->_write_dataset_data( $fg_dataset );
  $fg_dataset->data->process_metadata();

  return $fg_dataset;
}

sub get_bg_dataset {
  my $self = shift;
  my ( $fg_dataset ) = @_;

  my $root_dataset_id = $self->{query}->param('current_root_dataset_id');

#  my @pnames = $self->{query}->param;
#  my @dataset_ids = map { /chk_dataset_(\d+)/ } grep { /^chk_dataset/ } @pnames;
#  my @action_ids = map { /chk_action_(\d+)/ } grep { /^chk_action/ } @pnames;

#  if ( scalar @dataset_ids == 1 ) {
 #   $bg_dataset = EP::Common::DBI::Dataset->retrieve ( $dataset_ids[0] );
 # } elsif ( scalar @action_ids == 1 ) {
 #   $bg_dataset = EP::Common::DBI::Action->retrieve ( $action_ids[0] )->dst_dataset;
 # }

  # throw error if no background dataset is retrieved
  my $bg_dataset = EP::Common::DBI::Dataset->retrieve ( $root_dataset_id );
  throw EP::Common::Error ( -text => "Unable to retrieve background dataset" ) if !$bg_dataset;

  return $bg_dataset if $bg_dataset->type eq "seq";

  # Get the gene names ... for now just assume the NAMES entry in the row annotation is ok...
  my @ids = map { $bg_dataset->data->rows->{$_}{full_annot}->Values(0) }  sort { $a <=> $b } keys %{ $bg_dataset->data->rows } if $bg_dataset;
  # get the background sequences
  $self->{session}->status( "Retrieving background sequences ..." );

  # @ids - can try to get background sequences via another dataset...
#  my $bg_seq = EP::Core::Data::Parser::seq->_get_sequences ( $fg_dataset->species, \@ids, $self->{query}->param('upstream_bp') );  
  my $bg_seq = EP::Core::Data::Parser::seq->_get_sequences ( $fg_dataset->species, "", $self->{query}->param('upstream_bp') );
  $bg_dataset = $fg_dataset->copy ( { name => undef } );

  my $data = new EP::Core::Dataset::Data::seq ( sequences => $bg_seq );
  $bg_dataset->data( $data );
  $bg_dataset->data->start_bp( -1 * $self->{query}->param('upstream_bp') );
  $bg_dataset->data->end_bp( -1 );

  EP::Core::Data::Parser::seq->_write_dataset_data( $bg_dataset );
  $bg_dataset->data->process_metadata();

  return $bg_dataset;
}

sub _run_SPEXS {
  my $self = shift;

  throw EP::Common::Error ( -value => AFFY_TRANSFORM_RAW_DATA ) if $self->{src_dataset}->type eq "affy";

  my $fg_dataset = $self->get_fg_dataset();
  $self->{dst_dataset} = $fg_dataset;
  my $bg_dataset = $self->get_bg_dataset($fg_dataset);

  my (@params);

  my $dataset = $self->{src_dataset};
  my $folder = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $dataset->folder->location;

  # Minimum number of sequences that should have the pattern
  my $min_seq = int ( ( $self->{query}->param('-ms') / 100 ) * $fg_dataset->data->{nr_sequences} );

  # just in case something went wrong... set minimum number to 1 otherwise spexs will never finish
  $min_seq = 1 if $min_seq < 1;

  # parameter for foreground sequences
  my $fg_dataset_file = $folder . "/" . $fg_dataset->filename . ".fasta";
  push @params, "-f $fg_dataset_file -ms $min_seq";

  # parameter for background sequences
  my $bg_dataset_file = $folder . "/" . $bg_dataset->filename . ".fasta";
  push @params, "-f $bg_dataset_file -ms $min_seq";
  warn "Retrieved background sequences\n" if DEBUG;

  my $dataset_folder_path = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $dataset->folder->location;
  my $spexs_stdout = new File::Temp ( DIR => $dataset_folder_path,
				      UNLINK => 0,
				      SUFFIX => '.spexs' );

  # setup the rest of the parameters
  push @params, ("-showratio 1", "-freq_set 1", "-genorder 2");
  # SPEXS needs a config file to write temporary stuff to
  push @params, "-cf " . $EP::Config::h->{TMP} . "/spexs.$$.cf";
  foreach my $param ( qw( -goodratio -binomial_prob  ) ) {
    push @params, $param . " " . $self->{query}->param( $param );
  }

  my $jobs = { 0 => { program => $EP::Config::h->{SPEXS_CMD},
		      param => join( " ", @params ),
		      stdout => $spexs_stdout,
		      stderr => $spexs_stdout,
		    } };

  $self->{session}->status( "Running SPEXS ..." );
  die "Failed to run SPEXS!" if ( EP::Common::General::submit_jobs_and_wait ( $jobs ) ) != 1;
  my $results = $self->_visualize_SPEXS ( $spexs_stdout );

  return $results;
}

sub _visualize_SPEXS {
  my $self = shift;
  my $file = shift;

  my $patterns = $self->_read_SPEXS_patterns( $file );
  my $spexs_output = new File::Temp ( DIR => "$EP::Config::h->{EPNG_USRDATA_PATH}/" . $self->{src_dataset}->folder->location,
				      UNLINK => 0,
				      SUFFIX => ".spexs.xml" );

  print $spexs_output '<?xml version="1.0"?>';
  print $spexs_output '<table>';
  print $spexs_output '<head>';
  print $spexs_output '<data>Pattern</data>';
  print $spexs_output "<data>Cluster Frequency</data>";
  print $spexs_output "<data>Genome Frequency</data>";
  print $spexs_output "<data>Ratio</data>";
  print $spexs_output "<data>Binomial Probability</data>";
  print $spexs_output '</head>';

  # For now just sort on binomial probability, but could be ratio as well
  foreach my $pattern ( sort { $patterns->{$a}{binomial_prob} <=> $patterns->{$b}{binomial_prob} } keys %$patterns ) {
    print $spexs_output "<row>";
    print $spexs_output "<data>$_</data>" foreach $pattern, $patterns->{$pattern}{cluster_freq}, $patterns->{$pattern}{genome_freq},
			       $patterns->{$pattern}{ratio}, $patterns->{$pattern}{binomial_prob};
    print $spexs_output "</row>";
  }

  print $spexs_output '</table>';
  close $spexs_output;

  my $spexs_outfile = basename ( $spexs_output );

  return {  type  => "xml/table",
	    value => $EP::Config::h->{EPNG_USRDATA_URL} . "/" . $self->{src_dataset}->folder->location . "/" . $spexs_outfile,
	    label => "Motifs found by SPEXS" };
}

sub _read_SPEXS_patterns {
  my $self = shift;
  my $file = shift;
  local(*SPEXS);
  my ($patterns);

  open(SPEXS, $file) or die "Couldn't open $file, $!\n";
  while (<SPEXS>) {
    next if /^#|^\s/;
    chomp;
    if ( $_ =~ /^([ACGT]+)\s/ ) {
      # Found an actual pattern
      my $pattern = $1;
      my ($cluster_freq) = $_ =~ /1:(.+?)\s/;
      my ($genome_freq) = $_ =~ /2:(.+?)\s/;
      my ($ratio) = $_ =~ /R:([^\s]+)/;
      my ($prob) = $_ =~ /BP:([^\s]+)/;

      $patterns->{$pattern} = { cluster_freq => $cluster_freq,
				genome_freq => $genome_freq,
				ratio => $ratio,
				binomial_prob => $prob,
			      };
    }
  }
  close SPEXS;
    
  return $patterns;
}

1;
